
%%% The Ambiguity of Fishing for Fun, October 2022
%This function comes from "Econometric Analysis of Stochastic Dominance:
%Concepts, Methods, Tools, and Applications" (2001) by Y-Y Whang (p.
%223-224).

%--------------------------------------------------------Input-----------------------------------------------------------------%
function [p_value]=TestFSD(sample1,sample2,bmethod)

sample1=randn(200,1);
sample2=randn(200,1);
bmethod='Multiplier';
B=1000;
SDorder=1;
ngrid=100;

%---------------------------------------------------Subriutine-----------------------------------------------------------------%

%Measure the size of data
N1=size(sample1,1);
N2=size(sample2,1);

%Construct a support
pooled=sort([sample1;sample2]);
if ngrid==0
    grid=pooled;
else
    grid=linspace(min(pooled),max(pooled),ngrid);
end

%Compute the Test Statistic
operator=@(X,z)(X<=z).*(z-X).^(SDorder-1)/factorial(SDorder-1);
rawcdf1=mean(bsxfun(operator,sample1,reshape(grid, [1,1,ngrid])),1);
rawcdf2=mean(bsxfun(operator,sample2,reshape(grid, [1,1,ngrid])),1);
cdf1=squeeze(rawcdf1); %ECDF
cdf2=squeeze(rawcdf2); %ECDF
stat=sqrt(N1*N2/(N1+N2))*max(cdf1-cdf2); % test statistic

%Multiplier method
if strcmp(bmethod,'Multiplier')==1
    temp1=bsxfun(operator,sample1,reshape(grid,[1,1,ngrid]))-repmat(rawcdf1,[N1,1,1]);
    temp2=bsxfun(operator,sample2,reshape(grid,[1,1,ngrid]))-repmat(rawcdf2,[N2,1,1]);
    bcdf1=sqrt(N1)*mean(repmat(temp1,[1,B,1]).*randn(N1,B,ngrid),1);
    bcdf2=sqrt(N2)*mean(repmat(temp2,[1,B,1]).*randn(N2,B,ngrid),1);
    lambda=N2/(N1+N2);
    bksstat=max(sqrt(lambda)*bcdf1-sqrt(1-lambda)*bcdf2,[],3);
    
%Recentered bootstrap
elseif strcmp(bmethod,'Recentered')==1
    index1=randi(N1,N1,B);
    index2=randi(N2,N2,B);
    bsample1=sample1(index1); %bootstrap sample
    bsample2=sample2(index2); %bootstrap sample
    bcdf1=mean(bsxfun(operator,bsample1,reshape(grid,[1,1,ngrid])),1)-repmat(rawcdf1,[1,B,1]);
    bcdf2=mean(bsxfun(operator,bsample2,reshape(grid,[1,1,ngrid])),1)-repmat(rawcdf2,[1,B,1]);
    bksstat=sqrt(N1*N2/(N1+N2))*max(bcdf1-bcdf2,[],3);
    
%Pooled sample bootsrap
elseif strcmp(bmethod,'Pooled')==1
    index=randi(N2+N1,B);  %N1+N2<B
    index1=index(1:N1,:);
    index2=index(N1+1:N2+N1,:);
    bsample1=pooled(index1); %bootsrap sample
    bsample2=pooled(index2); %bootsrap sample
    bcdf1=mean(bsxfun(operator,bsample1,reshape(grid, [1,1,ngrid])),1);
    bcdf2=mean(bsxfun(operator,bsample2,reshape(grid, [1,1,ngrid])),1);
    bksstat=sqrt(N1*N2/(N1+N2))*max(bcdf1-bcdf2,[],3);
end
%--------------------------------------------Output----------------------------------------------------------------------------%

p_value=sum(bksstat>stat)/B;

end
    
    
    
    
    
    
    
    
